##### Dan Hopkins, 7/8/2008
###################################
###################################

### this code produces the analyses in 
### Hopkins, Daniel J.  "The Limited Impacts of Ethnic and Racial Diversity"
### Forthcoming, American Politics Research

### Please note that the code below allows estimation of results for both cities and counties.  
### The article and data only deal with U.S. cities.  

setwd("C:/Users/Dan/Documents/cityspnd")
load("cityspending031410.Rdata")
#load("citydata070808.Rdata")
#load("citydata052808.Rdata")

### 
CY2 <- CY12

### libarires
library(Matching)
library(MatchIt)
library(xtable)
library(Matrix)
library(MASS)
library(lme4)
library(arm)
library(mice)

### remove cities with no spending in one of the years in question
CY2 <- CY2[! (CY2$spendtot92==0 | CY2$spendtot==0),]

####
#### 1950 - 1960

CY2$PCNHBL50B <- CY2$PCNHBL50/100
CY2$MEDFAMINCD50S <- CY2$MEDFAMINC59-CY2$MEDFAMINC49
CY2$PCNHBLD50S <- CY2$PCNHBL60B-CY2$PCNHBL50B
CY2$OVER65D50S <- (CY2$OVER6560-CY2$OVER6550)/100
CY2$OVER65D60S <- CY2$OVER6570-(CY2$OVER6560)/100

CY2$crimeopsh73[CY2$crimeopsh73==0] <- NA
CY2$crimeopsh73[CY2$AREANAME.x=="Pico Rivera city"] <- CY2$crimesh73[CY2$AREANAME.x=="Pico Rivera city"] <- NA
### another potential outlier = Florrissant city"


#### levels

lout2<-lmer(crimesh50~PCNHBL50B+MEDFAMINC49+OVER6550+LGPOP50+(1|STATE),data=CY2)
lout2<-lmer(crimesh60~PCNHBL60B+MEDFAMINC59+OVER6560+LGPOP60+(1|STATE),data=CY2)
 
vars <- c("crimesh50","crimesh60","crimesh65","crimeopsh73","healthsh65","healthsh60","healthsh50","roadsh60","roadsh65","roadsh48","PCNHBL50B","PCNHBLD50S","PCNHBL70.x","MEANINC70","MEDFAMINCD50S","OVER6570","MEDFAMINC49","OVER6550","OVER65D50S","LGPOP50","LGPOPD50S","LGPOP70","PCNHWH70.x","PCPOOR70")
which(! vars %in% colnames(CY2))

XX<-subset(CY2,select=c("crimesh50","crimesh60","crimesh65","crimeopsh73","healthsh65","healthsh60","healthsh50","roadsh60","roadsh65","roadsh48","PCNHBL50B","PCNHBLD50S","PCNHBL70.x","MEANINC70","MEDFAMINCD50S","OVER6570","MEDFAMINC49","OVER6550","OVER65D50S","LGPOP50","LGPOPD50S","LGPOP70","PCNHWH70.x","PCPOOR70"))

X2 <- as.data.frame(XX)
mm <- mice(X2)

resmath <-resmatr <- resmatc <- matrix(NA,6,4)
colnames(resmath) <- colnames(resmatc) <- colnames(resmatr) <- c("beta.l","se.l","beta.d","se.d")

for(i in 1:5){
  CY3 <- complete(mm,action=i)
  CY3$STATE <- CY2$STATE 
  lout2<-lmer(crimesh60-crimesh50~PCNHBL50B+PCNHBLD50S+MEDFAMINC49+MEDFAMINCD50S+OVER6550+OVER65D50S+LGPOP50+LGPOPD50S+(1|STATE),data=CY3)
  resmatc[i,1] <- fixef(lout2)[2]
  resmatc[i,2] <- se.coef(lout2)$fixef[2]
  resmatc[i,3] <- fixef(lout2)[3]
  resmatc[i,4] <- se.coef(lout2)$fixef[3]

  lout1<-lmer(roadsh60-roadsh48~PCNHBL50B+PCNHBLD50S+MEDFAMINC49+MEDFAMINCD50S+OVER6550+OVER65D50S+LGPOP50+LGPOPD50S+(1|STATE),data=CY3)
  resmatr[i,1] <- fixef(lout1)[2]
  resmatr[i,2] <- se.coef(lout1)$fixef[2]
  resmatr[i,3] <- fixef(lout1)[3]
  resmatr[i,4] <- se.coef(lout1)$fixef[3]

  lout3 <- lmer(healthsh60-healthsh50~PCNHBL50B+PCNHBLD50S+MEDFAMINC49+MEDFAMINCD50S+OVER6550+OVER65D50S+LGPOP50+LGPOPD50S+(1|STATE),data=CY3)
  resmath[i,1] <- fixef(lout3)[2]
  resmath[i,2] <- se.coef(lout3)$fixef[2]
  resmath[i,3] <- fixef(lout3)[3]
  resmath[i,4] <- se.coef(lout3)$fixef[3]

}
resmatc60 <- resmatc

CYmat60 <- matrix(NA,6,2)
                  
source("C:/Users/Dan/Documents/Methods/miscRcode.R")
###top of CYmat73 = levels, crime then roads then healthcare
CYmat60[1,1] <- mean(resmatc[1:5,1])-1.96*mi.se(resmatc[1:5,1],resmatc[1:5,2])
CYmat60[1,2] <- mean(resmatc[1:5,1])+1.96*mi.se(resmatc[1:5,1],resmatc[1:5,2])

CYmat60[2,1] <- mean(resmatr[1:5,1])-1.96*mi.se(resmatr[1:5,1],resmatr[1:5,2])
CYmat60[2,2] <-  mean(resmatr[1:5,1])+1.96*mi.se(resmatr[1:5,1],resmatr[1:5,2])

CYmat60[3,1] <- mean(resmath[1:5,1])-1.96*mi.se(resmath[1:5,1],resmath[1:5,2])
CYmat60[3,2] <-  mean(resmath[1:5,1])+1.96*mi.se(resmath[1:5,1],resmath[1:5,2])

###changes for crime
CYmat60[4,1] <- mean(resmatc[1:5,3])-1.96*mi.se(resmatc[1:5,3],resmatc[1:5,4])
CYmat60[4,2] <-  mean(resmatc[1:5,3])+1.96*mi.se(resmatc[1:5,3],resmatc[1:5,4])

###changes for roads
CYmat60[5,1] <- mean(resmatr[1:5,3])-1.96*mi.se(resmatr[1:5,3],resmatr[1:5,4])
CYmat60[5,2] <-  mean(resmatr[1:5,3])+1.96*mi.se(resmatr[1:5,3],resmatr[1:5,4])

CYmat60[6,1] <- mean(resmath[1:5,3])-1.96*mi.se(resmath[1:5,3],resmath[1:5,4])
CYmat60[6,2] <-  mean(resmath[1:5,3])+1.96*mi.se(resmath[1:5,3],resmath[1:5,4])

rownames(CYmat60)  <- c("lev.cr","lev.rd","lev.he","ch.cr","ch.rd","ch.he")

CYmat60ad <- CYmat60
upper <- .9
lower <- .1

CYmat60ad[1:3,] <- CYmat60[1:3,]*(quantile(CY2$PCNHBL50B,upper,na.rm=T)- quantile(CY2$PCNHBL50B,lower,na.rm=T)) 
CYmat60ad[4:6,] <- CYmat60[4:6,]*(quantile(CY2$PCNHBLD50S,upper,na.rm=T)- quantile(CY2$PCNHBLD50S,lower,na.rm=T)) 

########################

XX<-subset(CY2,select=c("healthopsh73","healthsh60","healthsh65","roadsh60","roadopsh73","roadsh65","crimesh65","crimeopsh73","crimesh60","PCNHBL50B","PCNHBLD50S","PCNHBL70.x","MEANINC70","MEDFAMINCD50S","MEDFAMINC49","OVER6570","OVER6550","OVER65D50S","LGPOP50","LGPOP70","LGPOPD50S","PCNHWH70.x","PCPOOR70"))

X2 <- as.data.frame(XX)
mm <- mice(X2)

resmath <-resmatr <- resmatc <- matrix(NA,6,4)
colnames(resmath) <- colnames(resmatc) <- colnames(resmatr) <- c("beta.l","se.l","beta.d","se.d")

for(i in 1:5){
  CY3 <- complete(mm,action=i)
  CY3$STATE <- CY2$STATE
  lout2<-lmer(crimesh65-crimesh60~PCNHBL50B+PCNHBLD50S+MEDFAMINC49+MEDFAMINCD50S+OVER6550+OVER65D50S+LGPOP50+LGPOPD50S+(1|STATE),data=CY3)
  resmatc[i,1] <- fixef(lout2)[2]
  resmatc[i,2] <- se.coef(lout2)$fixef[2]
  resmatc[i,3] <- fixef(lout2)[3]
  resmatc[i,4] <- se.coef(lout2)$fixef[3]

  lout1<-lmer(roadsh65-roadsh60~PCNHBL50B+PCNHBLD50S+MEDFAMINC49+MEDFAMINCD50S+OVER6550+OVER65D50S+LGPOP50+LGPOPD50S+(1|STATE),data=CY3)
  resmatr[i,1] <- fixef(lout1)[2]
  resmatr[i,2] <- se.coef(lout1)$fixef[2]
  resmatr[i,3] <- fixef(lout1)[3]
  resmatr[i,4] <- se.coef(lout1)$fixef[3]

  lout3 <- lmer(healthsh65-healthsh60~PCNHBL50B+PCNHBLD50S+MEDFAMINC49+MEDFAMINCD50S+OVER6550+OVER65D50S+LGPOP50+LGPOPD50S+(1|STATE),data=CY3)
  resmath[i,1] <- fixef(lout3)[2]
  resmath[i,2] <- se.coef(lout3)$fixef[2]
  resmath[i,3] <- fixef(lout3)[3]
  resmath[i,4] <- se.coef(lout3)$fixef[3]
}
resmatc65 <- resmatc
CYmat65 <- matrix(NA,6,2)
                  
###top of CYmat73 = levels, crime then roads then healthcare
CYmat65[1,1] <- mean(resmatc[1:5,1])-1.96*mi.se(resmatc[1:5,1],resmatc[1:5,2])
CYmat65[1,2] <- mean(resmatc[1:5,1])+1.96*mi.se(resmatc[1:5,1],resmatc[1:5,2])

CYmat65[2,1] <- mean(resmatr[1:5,1])-1.96*mi.se(resmatr[1:5,1],resmatr[1:5,2])
CYmat65[2,2] <-  mean(resmatr[1:5,1])+1.96*mi.se(resmatr[1:5,1],resmatr[1:5,2])

CYmat65[3,1] <- mean(resmath[1:5,1])-1.96*mi.se(resmath[1:5,1],resmath[1:5,2])
CYmat65[3,2] <-  mean(resmath[1:5,1])+1.96*mi.se(resmath[1:5,1],resmath[1:5,2])

###changes for crime
CYmat65[4,1] <- mean(resmatc[1:5,3])-1.96*mi.se(resmatc[1:5,3],resmatc[1:5,4])
CYmat65[4,2] <-  mean(resmatc[1:5,3])+1.96*mi.se(resmatc[1:5,3],resmatc[1:5,4])

###changes for roads
CYmat65[5,1] <- mean(resmatr[1:5,3])-1.96*mi.se(resmatr[1:5,3],resmatr[1:5,4])
CYmat65[5,2] <-  mean(resmatr[1:5,3])+1.96*mi.se(resmatr[1:5,3],resmatr[1:5,4])

CYmat65[6,1] <- mean(resmath[1:5,3])-1.96*mi.se(resmath[1:5,3],resmath[1:5,4])
CYmat65[6,2] <-  mean(resmath[1:5,3])+1.96*mi.se(resmath[1:5,3],resmath[1:5,4])

rownames(CYmat65)  <- c("lev.cr","lev.rd","lev.he","ch.cr","ch.rd","ch.he")

CYmat65ad <- CYmat65
upper <- .9
lower <- .1

CYmat65ad[1:3,]<-CYmat65[1:3,]*(quantile(CY2$PCNHBL50B,upper,na.rm=T)- quantile(CY2$PCNHBL50B,lower,na.rm=T)) 
CYmat65ad[4:6,]<-CYmat65[4:6,]*(quantile(CY2$PCNHBLD50S,upper,na.rm=T)- quantile(CY2$PCNHBLD50S,lower,na.rm=T)) 


###
### plots--just crime

txt1 <- "Impact of Baseline \n Diversity: Cities"
txt2 <- "Shift from Homogeneous to Diverse Environment"
txt3 <- "Change in Spending Share"

postscript("C:/Users/Dan/Documents/cityspnd/spend60sdec09.eps",width=16)
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/graphs/crime60sjuly08.eps",width=16)
#pdf("/nfs/fs1/projects/dhopkins/cityspnd/graphs/f50t65regF.pdf",width=16)

par(mfcol=c(1,2))
rnn <- c("Criminal Justice","Roads","Health")
n <- 3
k <- 1.25
dist <- .06

mn <- min(c(CYmat60ad,CYmat65ad))-.01
mx <- max(c(CYmat60ad,CYmat65ad))+.01

#mn <- -.035
#mx <- 0.035

plot(1,type="n",xlim=c(mn,mx),ylim=c(1,2),main=txt1,sub=txt2,xlab=txt3,ylab="",yaxt="n",cex.lab=1.6,cex.main=1.8)
abline(v=0,col="grey")
for(i in c(1,2,3)){
  lines(x=c(CYmat60ad[i,1],CYmat60ad[i,2]),y=c(k,k),lwd=2)
  lines(x=c(CYmat65ad[i,1],CYmat65ad[i,2]),y=c(k-dist,k-dist),lty=5,lwd=2)

  points(x=mean(c(CYmat60ad[i,1],CYmat60ad[i,2])),y=k,pch=4,cex=2)
  text(rnn[i],x=mean(c(CYmat60ad[i,1],CYmat60ad[i,2])),y=k+dist*1.2,cex=2)

  points(x=mean(c(CYmat65ad[i,1],CYmat65ad[i,2])),y=k-dist,pch=4,cex=2)
#  text("1960-1965",x=mean(c(CYmat65ad[i,1],CYmat60ad[i,2])),y=k-dist+dist/4,cex=2)

  k <- k+.3
}

#initialize counter
k <- 1.25
txt1 <- "Impact of Increasing \n Diversity: Cities"
txt2 <- "Shift from Static to Rapidly Diversifying Environment"
txt3 <- "Change in Spending Share"


mn <- min(c(CYmat60ad,CYmat65ad))
mx <- max(c(CYmat60ad,CYmat65ad))

#mn <- -.035
#mx <- 0.035

plot(1,type="n",xlim=c(mn,mx),ylim=c(1,2),main=txt1,sub=txt2,xlab=txt3,ylab="",yaxt="n",cex.lab=1.6,cex.main=1.8)
abline(v=0,col="grey")
for(i in c(4,5,6)){
  lines(x=c(CYmat60ad[i,1],CYmat60ad[i,2]),y=c(k,k),lwd=2)
  lines(x=c(CYmat65ad[i,1],CYmat65ad[i,2]),y=c(k-dist,k-dist),lty=5,lwd=2)

  points(x=mean(c(CYmat60ad[i,1],CYmat60ad[i,2])),y=k,pch=4,cex=2)
  text(rnn[i-3],x=mean(c(CYmat60ad[i,1],CYmat60ad[i,2])),y=k+dist*1.2,cex=2)

  points(x=mean(c(CYmat65ad[i,1],CYmat65ad[i,2])),y=k-dist,pch=4,cex=2)
#  text("1960-1965",x=mean(c(CYmat65ad[i,1],CYmat60ad[i,2])),y=k-dist+dist/4,cex=2)

  k <- k+.3
}

dev.off()

#### Later resuls, data

load("C:/Users/Dan/Documents/cityspnd/county070808.Rdata")
tt6 <- tt7

####### Counties 1960s 
#######

tt10 <- subset(tt7,select=c("roadsh58","roadsh73","crimesh58","crimesh73","PCINSCH60","PCINSCH50","HERFD50S","HERF50","HERF70","COPCNHBL50","POP50","POP70","PCNHBL70"))
mm <- mice(tt10)

out<-lm.mids(crimesh73-crimesh58 ~ PCINSCH60+PCINSCH50+HERF50+HERFD50S,data=mm)
out<-lm.mids(roadsh73-roadsh58 ~ PCINSCH60+PCINSCH50+HERF50+HERFD50S,data=mm)

tt10 <- subset(tt7,select=c("roadsh83","roadsh73","crimesh83","crimesh73","PCINSCH60","HERF70","HERFD70S","HERF60","COPCNHBL60","POP60","POP70","PCNHBL80","MEANINC70","PCPOOR70","LGPOP70"))
library(mice)
mm <- mice(tt10)

out<-lm.mids(crimesh83-crimesh73 ~ HERF70+HERFD70S+POP70+PCPOOR70+MEANINC70,data=mm)
out<-lm.mids(roadsh83-roadsh73 ~ HERF70+HERFD70S,data=mm)

####MATCHING 1960
tt6$rownum <- 1:dim(tt6)[1]
dtt <- subset(tt6,select=c("HERFD60S","PCINSCH60","LGPOP60","BRTHST60","PCNHBL60","PCNHAS60","PCNHNA60","rownum"))
dt2 <- na.omit(dtt)
              
#TRvec <- 1*(dt2$HERFD60S<quantile(dt2$HERFD60S,na.rm=T,.3))
#mout2 <- matchit(TRvec ~ LGPOP60 + BRTHST60+PCINSCH60+PCNHBL60+PCNHAS60+PCNHNA60,data=dt2)
#matcheddta <- match.data(mout2,group="all")
#tt7 <- tt6[matcheddta$rownum,]

##############SUMMARIZE DATA
library(xtable)
####Summarize Changes 

cn4 <- c("roadsh73","crimesh73","healthsh73","sewerssh73","parkssh73","firesh73","librarysh73","housingsh73","transitsh73")
cn2 <- c("roadTD07S","crimeTD07S","healthTD07S","sewersTD07S","parksTD07S","fireTD07S","libraryTD07S","housingTD07S","transitTD07S")
#cn3 <- c("roadTD90S","crimeTD90S","healthTD90S","sewersTD90S","parksTD90S","fireTD90S","libraryTD90S","housingTD90S","transitTD90S")

rn <- c("Roads","Criminal Justice","Health","Sanitation","Parks","Fire","Libraries","Housing","Transit")

mat <- matrix(NA,length(cn2),4)
rownames(mat) <- rn

for(i in 1:length(cn2)){
    txt <- paste("mat[i,1] <-mean(CY2$",cn4[i],",na.rm=T)",sep="")
    eval(parse(text=txt))
    txt2 <- paste("mat[i,2] <- sd(CY2$",cn4[i],",na.rm=T)",sep="")
    eval(parse(text=txt2))
    txt <- paste("mat[i,3] <-mean(CY2$",cn2[i],",na.rm=T)",sep="")
    eval(parse(text=txt)) 
    txt2 <- paste("mat[i,4] <- sd(CY2$",cn2[i],",na.rm=T)",sep="")
    eval(parse(text=txt2))
}

colnames(mat) <-c("Mean, 73","SD","Mean, 73-02","SD")
rownames(mat) <- rn
xtable(mat,digits=rep(3,5))

###############how much spent
CY <- CY2
att2 <- c(mean(CY$transitsh,na.rm=T),mean(CY$roadsh,na.rm=T),mean(CY$healthsh,na.rm=T),mean(CY$crimesh,na.rm=T),mean(CY$parkssh,na.rm=T),mean(CY$housingsh,na.rm=T),mean(CY$firesh,na.rm=T),mean(CY$librarysh,na.rm=T),mean(CY$sewerssh,na.rm=T))
sum(att2)

att2 <- c(mean(tt6$transitsh,na.rm=T),mean(tt6$roadsh,na.rm=T),mean(tt6$healthsh,na.rm=T),mean(tt6$crimesh,na.rm=T),mean(tt6$parkssh,na.rm=T),mean(tt6$housingsh,na.rm=T),mean(tt6$firesh,na.rm=T),mean(tt6$librarysh,na.rm=T),mean(tt6$sewerssh,na.rm=T))
sum(att2)

###################Figure with spending shares

rn <- c("Transit","Roads","Healthcare","Criminal Justice","Parks","Housing","Fire","Libraries","Sanitation")

mattCO.zero <- c(sum(tt6$transitsh==0,na.rm=T),sum(tt6$roadsh==0,na.rm=T),sum(tt6$healthsh==0,na.rm=T),sum(tt6$crimesh==0,na.rm=T),sum(tt6$parkssh==0,na.rm=T),sum(tt6$housingsh==0,na.rm=T),sum(tt6$firesh==0,na.rm=T),sum(tt6$librarysh==0,na.rm=T),sum(tt6$sewerssh==0,na.rm=T))

mattCO <- c(mean(tt6$transitsh,na.rm=T),mean(tt6$roadsh,na.rm=T),mean(tt6$healthsh,na.rm=T),mean(tt6$crimesh,na.rm=T),mean(tt6$parkssh,na.rm=T),mean(tt6$housingsh,na.rm=T),mean(tt6$firesh,na.rm=T),mean(tt6$librarysh,na.rm=T),mean(tt6$sewerssh,na.rm=T))

names(mattCO) <- rn
matCO.sorted <-sort(mattCO,decreasing=T)

matt <- c(mean(CY$transitsh,na.rm=T),mean(CY$roadsh,na.rm=T),mean(CY$healthsh,na.rm=T),,mean(CY$crimesh,na.rm=T),mean(CY$parkssh,na.rm=T),mean(CY$housingsh,na.rm=T),mean(CY$firesh,na.rm=T),mean(CY$librarysh,na.rm=T),mean(CY$sewerssh,na.rm=T))

names(matt) <- rn
matt.sorted <- matt[order(mattCO,decreasing=T)]

png("/nfs/fs1/projects/dhopkins/cityspnd/spendingsum101106.png")
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/spendingsum91106.eps",height=9,width=9)
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/spendingsum.eps",height=6,width=6)
dotchart(matCO.sorted,main="Mean Share of City and County \n Spending by Category, 2002",col="red",xlim=c(0,.20),pch=15,ylab="",cex.main=1.8,cex.lab=1.8)
par(new=T)

dotchart(matt.sorted,main="",pch=8,xlim=c(0,.20),cex.lab=1.8)
legend(.15,8,c("County","City"),pch=c(15,8),col=c("red","black"),cex=1.2)
dev.off()

######SUMMARIZE KEY INDEPENDENT VARIABLES

###KEY CAUSAL VARIABLES
###COUNTY
ttc <- tt6[! (tt6$spendtot92 %in% c(0) | tt6$spendtot %in% c(0)),]

cn <-
  c("HERF00","HERFD07S","HERF70","LGPOP00","LGPOPD07S","LGPOP70","HIGINI00","HIGINID07S","HIGINI70f","PCPOOR00","PCPOORD07S","PCPOOR70","OVER6500","OVER65D07S","OVER6570","revigstpc","revigstpcD07S","revigstpc73","funcindex","dempct88","MEANINC00C","MEANINCD07S","MEANINC70C","MEHOVAL00")
    
rn <- c("Herf 00","$\Delta$ Herf 70-00","Herf
70","Lg. Pop. 00","$\Delta$ Lg. Pop. 70-00","Lg. Pop. 70","Inc. Gini
00","$\Delta$ Inc. Gini 70-00","Inc. Gini 70","Pct. Poor 00","$\Delta$
Pct. Poor 70-00","Pct. Poor 70","Pct. Over 65 00","$\Delta$ Pct. Over
65 70-00","Pct. Over 65 70","State Rev. per Cap. 02","$\Delta$ State
Rev. per Cap. 73-02","State Rev. per Cap. 73","Functional
Breadth","Pct. Dem. Vote 88","Mean Income 00","Mean Inc. 70-00","Mean
Inc. 70","Med. Home Val. 00")

mat <- matrix(NA,length(cn),6)
rownames(mat) <- rn

for(i in 1:length(cn)){
  txt <- paste("sout<-summary(ttc$",cn[i],")",sep="")
  eval(parse(text=txt))
  txt2 <- paste("sd.out <- sd(ttc$",cn[i],",na.rm=T)",sep="")
  eval(parse(text=txt2))
  mat[i,1] <-sout[1]
  mat[i,2] <-sout[3]
  mat[i,3] <-sout[4]
  mat[i,4] <-sout[6]
  mat[i,5] <- sd.out
  mat[i,6] <- sout[7]
  if(sout[7] %in% c(NA)) mat[i,6] <- 0
  
}

colnames(mat) <-c("Min","Median","Mean","Max","SD","Missing")
#write.table(mat,file="CYivsumm.csv",sep=",",na=".")

#mat2<-mat[order(mat[,3],decreasing=F),]
xtable(mat,digits=c(rep(2,(dim(mat)[2])),0))

############LEVELS

###KEY CAUSAL VARIABLES
###CITY

cn <-c("HERF00","HERFD07S","HERF70","LGPOP00","LGPOPD07S","LGPOP70","HIGINI00","HIGINID07S","HIGINI70f","PCPOOR00","PCPOORD07S","PCPOOR70","OVER6500","OVER65D07S","OVER6570","revigstpc","revigstpcD07S","revigstpc73","funcindex","WARD","CM","BDGTAUTH","MEANINC00C","MEANINCD07S","MEANINC70C","MEHOVAL00","crimerate91")
    
rn <- c("Herf 00","$\Delta$ Herf 70-00","Herf 70","Lg. Pop. 00","$\Delta$ Lg. Pop. 70-00","Lg. Pop. 70","Inc. Gini 00","$\Delta$ Inc. Gini 70-00","Inc. Gini 70","Pct. Poor
00","$\Delta$ Pct. Poor 70-00","Pct. Poor 70","Pct. Over 65 00","$\Delta$ Pct. Over 65 70-00","Pct. Over 65 70","State Rev. per
Cap. 02","$\Delta$ State Rev. per Cap. 73-02","State Rev. per Cap. 73","Functional Breadth","Ward
Elections","Council-Manager","Elected Budgted Authority","Mean Income 00","Mean Inc. 70-00","Mean Inc. 70","Med. Home
Val. 00","Crime Rate 1991")

mat <- matrix(NA,length(cn),6)
rownames(mat) <- rn

for(i in 1:length(cn)){
  txt <- paste("sout<-summary(CYC$",cn[i],")",sep="")
  eval(parse(text=txt))
  txt2 <- paste("sd.out <- sd(CYC$",cn[i],",na.rm=T)",sep="")
  eval(parse(text=txt2))
  mat[i,1] <-sout[1]
  mat[i,2] <-sout[3]
  mat[i,3] <-sout[4]
  mat[i,4] <-sout[6]
  mat[i,5] <- sd.out
  mat[i,6] <- sout[7]
  if(sout[7] %in% c(NA)) mat[i,6] <- 0
  
}

colnames(mat) <-c("Min","Median","Mean","Max","SD","Missing")
#write.table(mat,file="CYivsumm.csv",sep=",",na=".")

#mat2<-mat[order(mat[,3],decreasing=F),]
xtable(mat,digits=c(rep(2,(dim(mat)[2])),0))


#####DIVERSITY'S IMPACT ON EACH SPENDING AREA
###COUNTY 2002

library(lme4)
library(arm)

cn <-c("transitsh","roadsh","healthsh","crimesh","parkssh","housingsh","firesh","librarysh","sewerssh")

rn <- c("Transit","Roads","Healthcare","Crime","Parks","Housing","Fire","Libraries","Sewers")

mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn

for(i in 1:length(cn)){
  txt <- paste("lout<-lmer(",cn[i],"~HERF00+HIGINI00+MHINC00+MEHOVAL00+PCBCH00+PCPOOR00+OVER6500+PCURB90+LGPOP00+PCTOWN00+revigstpc+revigfedpc+funcindex+(1|STATE),data=tt6)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]
}

COMAT02 <- mat
#mat2<-mat[order(mat[,1]),]
#round(mat2,digits=4)

###COUNTY 1997
cn <- c("transitsh97","roadsh97","healthsh97","crimesh97","parkssh97","housingsh97","firesh97","librarysh97","sewerssh97")

mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn

for(i in 1:length(cn)){
  txt <- paste("lout<-lmer(",cn[i],"~HERF90+HIGINI90+MHINC90+MEHOVAL90+PCBCH90+PCPOOR90+OVER6590+PCURB90+LGPOP90+PCTOWN90+revigstpc+revigfedpc+funcindex+(1|STATE),data=tt6)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]

}

COMAT97 <- mat
mat2<-mat[order(mat[,1]),]
round(mat2,digits=4)


###COUNTY 1992
cn <- c("transitsh92","roadsh92","healthsh92","crimesh92","parkssh92","housingsh92","firesh92","librarysh92","sewerssh92")

mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn

for(i in 1:length(cn)){
  txt <- paste("lout<-lmer(",cn[i],"~HERF90+HIGINI90+MHINC90+MEHOVAL90+PCBCH90+PCPOOR90+OVER6590+PCURB90+LGPOP90+PCTOWN90+revigstpc+revigfedpc+funcindex+(1|STATE),data=tt6)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]
}

COMAT92 <- mat
mat2<-mat[order(mat[,1]),]
round(mat2,digits=4)


###COUNTY 1973
cn <- c("transitsh73","roadsh73","healthsh73","crimesh73","parkssh73","housingsh73","firesh73","librarysh73","sewerssh73")

mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn

for(i in 1:length(cn)){
  txt <-
    paste("lout<-lmer(",cn[i],"~HERF70+MEANINC70+HIGINI70f+OVER6570+PCPOOR70+LGPOP70+(1|STATE),data=tt6)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]
}

COMAT73 <- mat
mat2<-mat[order(mat[,1]),]
round(mat2,digits=4)

###############CITIES



################
CYC <- CY <- CY2 
CYC$HERFB00 <- CYC$PCNHBL00^2+CYC$PCHIS00^2+CYC$PCNHWH00^2#+CYC$PCNHNA00^2+CYC$PCNHAS00^2
CYC$HERFA00 <- CYC$PCNHBL00^2+CYC$PCHIS00^2+CYC$PCNHWH00^2+CYC$PCNHAS00^2#+CYC$PCNHNA00^2



cn <-c("transitsh","roadsh","healthsh","crimesh","parkssh","housingsh","firesh","librarysh","sewerssh")

rn <-c("Transit","Roads","Healthcare","Crime","Parks","Housing","Fire","Libraries","Sewers")

mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn

for(i in 1:length(cn)){
  txt <- paste("CYD <- CYC[CYC$",cn[i],">0,]",sep="")
  eval(parse(text=txt))

  txt <- paste("lout<-lmer(",cn[i],"~HERF00+HIGINI00+MHINC00+MEHOVAL00+PCBCH00+PCPOOR00+OVER6500+PCURB90+LGPOP00+PCTOWN00+revigstpc+revigfedpc+funcindex+(1|STATE),data=CYD)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]
}

CYMAT02 <- mat
#mat2<-mat[order(mat[,1]),]
#round(mat2,digits=4)

###CITY 1997
cn <-
  c("transitsh97","roadsh97","healthsh97","crimesh97","parkssh97","housingsh97","firesh97","librarysh97","sewerssh97")

mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn

for(i in 1:length(cn)){

  txt <- paste("CYD <- CYC[CYC$",cn[i],">0,]",sep="")
  eval(parse(text=txt))

  txt <- paste("lout<-lmer(",cn[i],"~HERF90+HIGINI90+MHINC90+MEHOVAL90+PCBCH90+PCPOOR90+OVER6590+PCURB90+LGPOP90+PCTOWN90+revigstpc+revigfedpc+funcindex+(1|STATE),data=CYD)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]

}

CYMAT97 <- mat
mat2<-mat[order(mat[,1]),]
round(mat2,digits=4)


###CITY 1992
cn <-
  c("transitsh92","roadsh92","healthsh92","crimesh92","parkssh92","housingsh92","firesh92","librarysh92","sewerssh92")



mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn

for(i in 1:length(cn)){

  txt <- paste("CYD <- CYC[CYC$",cn[i],">0,]",sep="")
  eval(parse(text=txt))
  txt <- paste("lout<-lmer(",cn[i],"~HERF90+HIGINI90+MHINC90+MEHOVAL90+PCBCH90+PCPOOR90+OVER6590+PCURB90+LGPOP90+PCTOWN90+revigstpc+revigfedpc+funcindex+(1|STATE),data=CYD)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]
}

CYMAT92 <- mat
mat2<-mat[order(mat[,1]),]
round(mat2,digits=4)


###CITY 1973
cn <-
  c("transitsh73","roadsh73","healthsh73","crimesh73","parkssh73","housingsh73","firesh73","librarysh73","sewerssh73")

mat <- matrix(NA,length(cn),2)
rownames(mat) <- rn


for(i in 1:length(cn)){
  txt <- paste("CYD <- CYC[CYC$",cn[i],">0,]",sep="")
  eval(parse(text=txt))
  txt <-
    paste("lout<-lmer(",cn[i],"~HERF70+MEANINC70+HIGINI70f+OVER6570+PCPOOR70+LGPOP70+(1|STATE),data=CYD)",sep="")
  eval(parse(text=txt))
  mat[i,1] <-fixef(lout)[2]
  mat[i,2] <-se.fixef(lout)[2]
}

CYMAT73 <- mat
mat2<-mat[order(mat[,1]),]
round(mat2,digits=4)

####print all results
#mat11<-round(cbind(CYMAT92,CYMAT97,CYMAT02),digits=3)
#colnames(mat11) <- c("$Beta$ Herf 92","SE","$Beta$ Herf 97","SE","$Beta$ Herf 02","SE")
#mat22<-round(cbind(COMAT92,COMAT97,COMAT02),digits=3)
#colnames(mat22) <- c("$Beta$ Herf 92","SE","$Beta$ Herf 97","SE","$Beta$ Herf 02","SE")
#xtable(mat11,digits=rep(3,1+dim(mat11)[2]))
#xtable(mat22,digits=rep(3,1+dim(mat11)[2]))

###graphical representation
mat11<-round(cbind(CYMAT73,CYMAT92,CYMAT97,CYMAT02),digits=6)
mat22<-round(cbind(COMAT73,COMAT92,COMAT97,COMAT02),digits=6)

mat11A <- mat11
#mat11A[10:11,] <- mat11[10:11,]/10

mat22A <- mat22
#mat22A[10:11,] <- mat22[10:11,]/10

mat33 <- matrix(NA,nrow(mat11),ncol(mat11))
mat33[,1] <- -mat11A[,1]-1.96*mat11A[,2]
mat33[,2] <- -mat11A[,1]+1.96*mat11A[,2]
mat33[,3] <- -mat11A[,3]-1.96*mat11A[,4]
mat33[,4] <- -mat11A[,3]+1.96*mat11A[,4]
mat33[,5] <- -mat11A[,5]-1.96*mat11A[,6]
mat33[,6] <- -mat11A[,5]+1.96*mat11A[,6]
mat33[,7] <- -mat11A[,7]-1.96*mat11A[,8]
mat33[,8] <- -mat11A[,7]+1.96*mat11A[,8]

mat44 <- matrix(NA,nrow(mat11),ncol(mat11))
mat44[,1] <- -mat22A[,1]-1.96*mat22A[,2]
mat44[,2] <- -mat22A[,1]+1.96*mat22A[,2]
mat44[,3] <- -mat22A[,3]-1.96*mat22A[,4]
mat44[,4] <- -mat22A[,3]+1.96*mat22A[,4]
mat44[,5] <- -mat22A[,5]-1.96*mat22A[,6]
mat44[,6] <- -mat22A[,5]+1.96*mat22A[,6]
mat44[,7] <- -mat22A[,7]-1.96*mat22A[,8]
mat44[,8] <- -mat22A[,7]+1.96*mat22A[,8]

n <- dim(mat44)[1]

rn <- c("Transit","Roads","Health","Crime","Parks","Housing","Fire","Libraries","Sanitation")
txttt <- "Increased Spending                  Reduced Spending"
#png("/nfs/fs1/projects/dhopkins/cityspnd/graphs/regCO10.png")
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/regCO80607.eps",height=8,width=8)
#pdf("/nfs/fs1/projects/dhopkins/cityspnd/regCO82006.pdf",height=8,width=8)
plot(y=c(1,dim(mat44)[1]*1.51),x=c(min(mat44,na.rm=T)-.01,max(mat44,na.rm=T)),type="n",yaxt="n",main=
     "Impact of Diversity on \n County Spending Shares by Year",ylab=txttt,xlab=
     "Conditional Impact of Level of Diversity",sub="(1992 Point Estimate in Parentheses)",cex.main=1.8,cex.lab=1.6)

mat23A <- round(mat22A,digits=3)

k <- 1
se <- order(mat22A[,3])
for(i in se){
  rn2 <- paste(rn[i]," (",-mat23A[i,3],")",sep="")
  lines(c(mat44[i,1],mat44[i,2]),c(k+.3,k+.3),lty=2)
  points(y=k+.3,x=mean(c(mat44[i,1],mat44[i,2])),pch=4,cex=.75)  
  if(min(mat44[i,],na.rm=T)>0 | i==3){
    text(x=min(mat44[i,],na.rm=T)-.045,y=k,rn2,cex=1.2)
  }else{
    text(x=max(mat44[i,],na.rm=T)+.065,y=k,rn2,cex=1.2)
  }
  lines(c(mat44[i,3],mat44[i,4]),c(k+.1,k+.1),lty=3)
  points(y=k+.1,x=mean(c(mat44[i,3],mat44[i,4])),pch=4,cex=.75)

  lines(c(mat44[i,5],mat44[i,6]),c(k-.1,k-.1),lty=1)
  points(y=k-.1,x=mean(c(mat44[i,5],mat44[i,6])),pch=4,cex=.75)
  lines(c(mat44[i,7],mat44[i,8]),c(k-.3,k-.3),lty=4)
  points(y=k-.3,x=mean(c(mat44[i,7],mat44[i,8])),pch=4,cex=.75)
  k <- k+1.5
}
legend(-.12,8,legend=c("1973","1992","1997","2002"),lty=c(2,3,1,4),cex=1.2)
abline(v=0,col="red")
dev.off()
#################
#################CITIES

mmn <- "Impact of Diversity on \n City Spending Shares by Year"
#png("/nfs/fs1/projects/dhopkins/cityspnd/graphs/regCY101106.png")

postscript("C:/Users/Dan/Documents/cityspnd/regCY071508.eps",height=8,width=8)
#pdf("/nfs/fs1/projects/dhopkins/cityspnd/graphs/regCY82006.pdf")
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/regCY052807.eps",height=8,width=8)
plot(y=c(1,dim(mat33)[1]*1.51),x=c(min(mat33,na.rm=T)-.01,max(mat33,na.rm=T)+.23),type="n",yaxt="n",main=mmn,ylab=txttt,xlab=
     "Conditional Impact of Level of Diversity",sub="(1992 Point Estimate in Parentheses)",cex.main=1.8,cex.lab=1.6)

mat13A <- round(mat11A,digits=3)

k <- 1
se <- order(mat11A[,3])
for(i in se){
  rn2 <- paste(rn[i]," (",-mat13A[i,3],")",sep="")
  lines(c(mat33[i,1],mat33[i,2]),c(k+.3,k+.3),lty=2)
  points(y=k+.3,x=mean(c(mat33[i,1],mat33[i,2])),pch=4,cex=.75)  
  
#if(min(mat33[i,],na.rm=T)>0){
#text(x=min(mat33[i,],na.rm=T)-.12,y=k,rn2,cex=1.2)
#}else{
  text(x=max(mat33[i,],na.rm=T)+.08,y=k,rn2,cex=1.2)
#}
  lines(c(mat33[i,3],mat33[i,4]),c(k+.1,k+.1),lty=3)
  points(y=k+.1,x=mean(c(mat33[i,3],mat33[i,4])),pch=4,cex=.75) 

  lines(c(mat33[i,5],mat33[i,6]),c(k-.1,k-.1),lty=1)
  points(y=k-.1,x=mean(c(mat33[i,5],mat33[i,6])),pch=4,cex=.75) 

  lines(c(mat33[i,7],mat33[i,8]),c(k-.3,k-.3),lty=4)
  points(y=k-.3,x=mean(c(mat33[i,7],mat33[i,8])),pch=4,cex=.75) 
  k <- k+1.5
}
legend(.25,6.2,legend=c("1973","1992","1997","2002"),lty=c(2,3,1,4),cex=1.2)
abline(v=0,col="red")
dev.off()

#################
################Make changes graph
################

CYmat <- matrix(NA,5,12)
colnames(CYmat) <- c("Lq02.90","Lq50.90","Lq97.90","Cq02.90","Cq50.90","Cq97.90","Lq02.97","Lq50.97","Lq97.97","Cq02.97","Cq50.97","Cq97.97")
cn <- rownames(CYmat) <- c("sewersTD90S","roadD90S","healthTD90S","housingTD90S","crimeTD90S")
cn2 <- c("sewersD97S","roadsD97S","healthTD97S","housingTD97S","crimeTD97S")
cn3 <- c("sewersD07S","roadsD07S","healthTD07S","housingTD07S","crimeTD07S")
#coefmat97 <- matrix(NA,5,length(cn))
#coefmat90 <- matrix(NA,5,length(cn))

clout1L<-lm(sewersTD90S~HERF90+HERFD90S+MHINC90+MHINCD90S+MEHOVAL90+MEHOVALD90S+OVER6590+OVER65D90S+LGPOP90+LGPOPD90S+PCPOOR90+PCPOORD90S+HIGINI90+HIGINID90S+PCTOWN90+PCTOWND90S+PCPA90+PCPAD90S+revigstpc92+revigstpcD90S,data=CYC)

mm1<-model.matrix(clout1L,data=CY)
mm2lc <- mm2ll <- mm2hc <- mm2hl <- mm2 <- apply(mm1,2,median)

mm2ll[2] <- quantile(CYC$HERF90,.2,na.rm=T)
mm2hl[2] <- quantile(CYC$HERF90,.8,na.rm=T)

mm2lc[3] <- quantile(CYC$HERFD90S,.2,na.rm=T)
mm2hc[3] <- quantile(CYC$HERFD90S,.8,na.rm=T)

##################MAKE 1970-1990 MATRIX

clout3L<-lm(sewersTD07S~HERF70+HERFD07S+MEANINC70+MEANINCD07S+OVER6570+OVER65D07S+LGPOP70+LGPOPD07S+PCPOOR70+PCPOORD07S+HIGINI70f+HIGINID07S+revigstpc73+revigstpcD07S,data=CYC)

m1<-model.matrix(clout3L,data=CY)
m2lc <- m2ll <- m2hc <- m2hl <- m2 <- apply(m1,2,median)
m2ll[2] <- quantile(CYC$HERF70,.2,na.rm=T)
m2hl[2] <- quantile(CYC$HERF70,.8,na.rm=T)

m2lc[3] <- quantile(CYC$HERFD07S,.2,na.rm=T)
m2hc[3] <- quantile(CYC$HERFD07S,.8,na.rm=T)

##############################################

CYmatA <- array(dim=c(5,12,5))

for(i in 1:length(cn)){
  txt <- paste("CYD <- CYC[CYC$",cn[i],">0,]",sep="")
  eval(parse(text=txt))  

### for above 10% black
#  txt <- paste("CYD <- CYC[CYC$",cn[i],">0 & CYC$PCNHBL00 > .1,]",sep="")
#  eval(parse(text=txt))  

  txt1 <- paste("clout1 <-lmer(",cn[i],"~HERF90+HERFD90S+MHINC90+MHINCD90S+MEHOVAL90+MEHOVALD90S+OVER6590+OVER65D90S+LGPOP90+LGPOPD90S+PCPOOR90+PCPOORD90S+HIGINI90+HIGINID90S+PCTOWN90+PCTOWND90S+PCPA90+PCPAD90S+revigstpc92+revigstpcD90S+(1|STATE),data=CYD,showCorrelation=F,hessian=T)",sep="")
  eval(parse(text=txt1))
  n <- 5000
  mu <- fixef(clout1)
  Sigma <- vcov(clout1)
  coefs <- mvrnorm(n,mu,Sigma)
  levels <-  coefs%*%mm2ll - coefs%*%mm2hl
  changes <- coefs%*%mm2lc - coefs%*%mm2hc
  CYmat[i,1:3] <-quantile(levels,c(0.025,.5,.975))
  CYmat[i,4:6] <- quantile(changes,c(0.025,.5,.975))
  
	XXX <- subset(CYD,select=c(cn3[i],"HERF70","HERFD07S","MEANINC70","MEANINCD07S","OVER6570","OVER65D07S","LGPOP70","LGPOPD07S","PCPOOR70","PCPOORD07S","HIGINI70f","HIGINID07S","revigstpc73","revigstpcD07S","STATE"))
#	txtt <- paste("XIV <- XXX[! XXX$",cn3[i],"==0,]",sep="")
#	eval(parse(text=txtt))
	MM <- mice(XXX)  
	holdveccha <- holdveclev <- c()
	for(j in 1:5){
		CYE <- complete(MM,action=j)
		txt1 <- paste("clout2<-lmer(",cn3[i],"~HERF70+HERFD07S+MEANINC70+MEANINCD07S+OVER6570+OVER65D07S+LGPOP70+LGPOPD07S+PCPOOR70+PCPOORD07S+HIGINI70f+HIGINID07S+revigstpc73+revigstpcD07S+(1|STATE),data=CYE,showCorrelation=F,hessian=T)",sep="")
  		eval(parse(text=txt1))
  		n <- 5000
  		mu <- fixef(clout2)
  		Sigma <- vcov(clout2)
  		coefs <- mvrnorm(n,mu,Sigma)
  		holdveclev <-  c(holdveclev,(coefs%*%m2ll - coefs%*%m2hl))
  		holdveccha <- c(holdveccha,(coefs%*%m2lc - coefs%*%m2hc))
	}
	CYmat[i,7:9] <-quantile(holdveclev,c(0.025,.5,.975))
  	CYmat[i,10:12] <- quantile(holdveccha,c(0.025,.5,.975))
}

round(CYmat[,1:6],digits=3)
round(CYmat[,7:12],digits=3)

####Appendix Table

i <- 5
txt1 <- paste("clout2<-lmer(",cn3[i],"~HERF70+HERFD07S+MEANINC70+MEANINCD07S+OVER6570+OVER65D07S+LGPOP70+LGPOPD07S+PCPOOR70+PCPOORD07S+HIGINI70f+HIGINID07S+revigstpc73+revigstpcD07S+(1|STATE),data=CY,showCorrelation=F,hessian=T)",sep="")
eval(parse(text=txt1))
display(clout2)

cln2 <- names(fixef(clout2))
mttt <- matrix(NA,length(cln2),2)
for(i in 1:15){
  mttt[i,1] <- fixef(clout2)[i]
  mttt[i,2] <- se.fixef(clout2)[i]
}
rownames(mttt) <- cln2
mttt2 <- mttt
mttt2[4:5,] <- mttt[4:5,]*1000000

xtable(mttt2,digits=c(0,3,3))


##################
####COmat2 <- COmat[,7:12]
CYmat2 <- CYmat[,7:12]
rnn <- c("Sanitation","Roads","Health","Housing","Criminal Justice")

#initialize counter
k <- 1
j <- 1.15
txt1 <- "Impact of Baseline \n Diversity: Cities"
txt2 <- "Shift from Homogeneous to Diverse Environment"
txt3 <- "Predicted Change in Spending Share"


postscript("C:/Users/Dan/Documents/cityspnd/sevenCI071508.eps",height=8,width=12)

#png("/nfs/fs1/projects/dhopkins/cityspnd/graphs/CYchanges101206.png",height=600,width=900)
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/graphs/sevenCI052808.eps",height=8,width=12)
#pdf("/nfs/fs1/projects/dhopkins/cityspnd/graphs/sevenCI.pdf",height=8,width=16)
par(mfcol=c(1,2))

n <- 5
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/graphs/levelCYff.eps",height=9,width=9)
plot(1,type="n",xlim=c(min(CYmat2),max(CYmat2)),ylim=c(1,n+.5),main=txt1,sub=txt2,xlab=txt3,ylab="",yaxt="n",cex.lab=1.6,cex.main=1.8)
abline(v=0,col="grey")
for(i in 1:n){
  lines(x=c(CYmat2[i,1],CYmat2[i,3]),y=c(k,k))
  points(x=CYmat2[i,2],y=k,pch=4)
  text(rnn[i],x=CYmat2[i,2],y=k+.3,cex=1.2)
  j <- j+1
  k <- k+1
}
#legend(legend=c("1973-1992","1992-2002"),x=-.029,y=4,lty=c(2,1),cex=1.2)
rnn <- c("Sanitation","Roads","Health","Housing","Criminal Justice")

#initialize counter
k <- 1
j <- 1.15
txt1 <- "Impact of Increasing \n Diversity: Cities"
txt2 <- "Shift from Slowly Diversifying to Rapidly Diversifying Environment"
txt3 <- "Predicted Change in Spending Share"

n <- 5
plot(1,type="n",xlim=c(min(CYmat2),max(CYmat2)),ylim=c(1,n+.5),main=txt1,sub=txt2,xlab=txt3,ylab="",yaxt="n",cex.lab=1.6,cex.main=1.8)
abline(v=0,col="grey")
for(i in 1:n){
  lines(x=c(CYmat2[i,4],CYmat2[i,6]),y=c(k,k))
  points(x=CYmat2[i,5],y=k,pch=4)
  text(rnn[i],x=CYmat2[i,5],y=k+.3,cex=1.2)
  j <- j+1
  k <- k+1
}

dev.off()

postscript("/nfs/fs1/projects/dhopkins/cityspnd/graphs/sevenCO80607.eps",height=8,width=12)
#pdf("/nfs/fs1/projects/dhopkins/cityspnd/graphs/seventiesCO.pdf",height=8,width=16)
par(mfcol=c(1,2))

########counties
rnn <- c("Sanitation","Roads","Health","Fire","Criminal Justice")
k <- 1
j <- 1.15
txt1 <- "Impact of Baseline \n Diversity: Counties"
txt2 <- "Shift from Homogeneous to Diverse Environment"
txt3 <- "Predicted Change in Spending Share"

n <- 5
#postscript("/nfs/fs1/projects/dhopkins/cityspnd/graphs/levelCYff.eps",height=9,width=9)
plot(1,type="n",xlim=c(min(COmat2),max(COmat2)),ylim=c(1,n+.5),main=txt1,sub=txt2,xlab=txt3,ylab="",yaxt="n",cex.lab=1.6,cex.main=1.8)
abline(v=0,col="grey")
for(i in 1:n){
  lines(x=c(COmat2[i,1],COmat2[i,3]),y=c(k,k))
  points(x=COmat2[i,2],y=k,pch=4)
  text(rnn[i],x=COmat2[i,2],y=k+.3,cex=1.2)
  j <- j+1
  k <- k+1
}
#legend(legend=c("1973-1992","1992-2002"),x=-.029,y=4,lty=c(2,1),cex=1.2)
rnn <- c("Sanitation","Roads","Health","Fire","Criminal Justice")

#initialize counter
k <- 1
j <- 1.15
txt1 <- "Impact of Increasing \n Diversity: Counties"
txt2 <- "Shift from Slowly Diversifying to Rapidly Diversifying Environment"
txt3 <- "Predicted Change in Spending Share"

n <- 5
plot(1,type="n",xlim=c(min(COmat2),max(COmat2)),ylim=c(1,n+.5),main=txt1,sub=txt2,xlab=txt3,ylab="",yaxt="n",cex.lab=1.6,cex.main=1.8)
abline(v=0,col="grey")
for(i in 1:n){
  lines(x=c(COmat2[i,4],COmat2[i,6]),y=c(k,k))
  points(x=COmat2[i,5],y=k,pch=4)
  text(rnn[i],x=COmat2[i,5],y=k+.3,cex=1.2)
  j <- j+1
  k <- k+1
}


dev.off()

